home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / lapack / dsteqr.f < prev    next >
Text File  |  1997-06-25  |  14KB  |  502 lines
  1.       SUBROUTINE DSTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO )
  2. *
  3. *  -- LAPACK routine (version 2.0) --
  4. *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
  5. *     Courant Institute, Argonne National Lab, and Rice University
  6. *     September 30, 1994
  7. *
  8. *     .. Scalar Arguments ..
  9.       CHARACTER          COMPZ
  10.       INTEGER            INFO, LDZ, N
  11. *     ..
  12. *     .. Array Arguments ..
  13.       DOUBLE PRECISION   D( * ), E( * ), WORK( * ), Z( LDZ, * )
  14. *     ..
  15. *
  16. *  Purpose
  17. *  =======
  18. *
  19. *  DSTEQR computes all eigenvalues and, optionally, eigenvectors of a
  20. *  symmetric tridiagonal matrix using the implicit QL or QR method.
  21. *  The eigenvectors of a full or band symmetric matrix can also be found
  22. *  if DSYTRD or DSPTRD or DSBTRD has been used to reduce this matrix to
  23. *  tridiagonal form.
  24. *
  25. *  Arguments
  26. *  =========
  27. *
  28. *  COMPZ   (input) CHARACTER*1
  29. *          = 'N':  Compute eigenvalues only.
  30. *          = 'V':  Compute eigenvalues and eigenvectors of the original
  31. *                  symmetric matrix.  On entry, Z must contain the
  32. *                  orthogonal matrix used to reduce the original matrix
  33. *                  to tridiagonal form.
  34. *          = 'I':  Compute eigenvalues and eigenvectors of the
  35. *                  tridiagonal matrix.  Z is initialized to the identity
  36. *                  matrix.
  37. *
  38. *  N       (input) INTEGER
  39. *          The order of the matrix.  N >= 0.
  40. *
  41. *  D       (input/output) DOUBLE PRECISION array, dimension (N)
  42. *          On entry, the diagonal elements of the tridiagonal matrix.
  43. *          On exit, if INFO = 0, the eigenvalues in ascending order.
  44. *
  45. *  E       (input/output) DOUBLE PRECISION array, dimension (N-1)
  46. *          On entry, the (n-1) subdiagonal elements of the tridiagonal
  47. *          matrix.
  48. *          On exit, E has been destroyed.
  49. *
  50. *  Z       (input/output) DOUBLE PRECISION array, dimension (LDZ, N)
  51. *          On entry, if  COMPZ = 'V', then Z contains the orthogonal
  52. *          matrix used in the reduction to tridiagonal form.
  53. *          On exit, if INFO = 0, then if  COMPZ = 'V', Z contains the
  54. *          orthonormal eigenvectors of the original symmetric matrix,
  55. *          and if COMPZ = 'I', Z contains the orthonormal eigenvectors
  56. *          of the symmetric tridiagonal matrix.
  57. *          If COMPZ = 'N', then Z is not referenced.
  58. *
  59. *  LDZ     (input) INTEGER
  60. *          The leading dimension of the array Z.  LDZ >= 1, and if
  61. *          eigenvectors are desired, then  LDZ >= max(1,N).
  62. *
  63. *  WORK    (workspace) DOUBLE PRECISION array, dimension (max(1,2*N-2))
  64. *          If COMPZ = 'N', then WORK is not referenced.
  65. *
  66. *  INFO    (output) INTEGER
  67. *          = 0:  successful exit
  68. *          < 0:  if INFO = -i, the i-th argument had an illegal value
  69. *          > 0:  the algorithm has failed to find all the eigenvalues in
  70. *                a total of 30*N iterations; if INFO = i, then i
  71. *                elements of E have not converged to zero; on exit, D
  72. *                and E contain the elements of a symmetric tridiagonal
  73. *                matrix which is orthogonally similar to the original
  74. *                matrix.
  75. *
  76. *  =====================================================================
  77. *
  78. *     .. Parameters ..
  79.       DOUBLE PRECISION   ZERO, ONE, TWO, THREE
  80.       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
  81.      $                   THREE = 3.0D0 )
  82.       INTEGER            MAXIT
  83.       PARAMETER          ( MAXIT = 30 )
  84. *     ..
  85. *     .. Local Scalars ..
  86.       INTEGER            I, ICOMPZ, II, ISCALE, J, JTOT, K, L, L1, LEND,
  87.      $                   LENDM1, LENDP1, LENDSV, LM1, LSV, M, MM, MM1,
  88.      $                   NM1, NMAXIT
  89.       DOUBLE PRECISION   ANORM, B, C, EPS, EPS2, F, G, P, R, RT1, RT2,
  90.      $                   S, SAFMAX, SAFMIN, SSFMAX, SSFMIN, TST
  91. *     ..
  92. *     .. External Functions ..
  93.       LOGICAL            LSAME
  94.       DOUBLE PRECISION   DLAMCH, DLANST, DLAPY2
  95.       EXTERNAL           LSAME, DLAMCH, DLANST, DLAPY2
  96. *     ..
  97. *     .. External Subroutines ..
  98.       EXTERNAL           DLAE2, DLAEV2, DLARTG, DLASCL, DLASET, DLASR,
  99.      $                   DLASRT, DSWAP, XERBLA
  100. *     ..
  101. *     .. Intrinsic Functions ..
  102.       INTRINSIC          ABS, MAX, SIGN, SQRT
  103. *     ..
  104. *     .. Executable Statements ..
  105. *
  106. *     Test the input parameters.
  107. *
  108.       INFO = 0
  109. *
  110.       IF( LSAME( COMPZ, 'N' ) ) THEN
  111.          ICOMPZ = 0
  112.       ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
  113.          ICOMPZ = 1
  114.       ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
  115.          ICOMPZ = 2
  116.       ELSE
  117.          ICOMPZ = -1
  118.       END IF
  119.       IF( ICOMPZ.LT.0 ) THEN
  120.          INFO = -1
  121.       ELSE IF( N.LT.0 ) THEN
  122.          INFO = -2
  123.       ELSE IF( ( LDZ.LT.1 ) .OR. ( ICOMPZ.GT.0 .AND. LDZ.LT.MAX( 1,
  124.      $         N ) ) ) THEN
  125.          INFO = -6
  126.       END IF
  127.       IF( INFO.NE.0 ) THEN
  128.          CALL XERBLA( 'DSTEQR', -INFO )
  129.          RETURN
  130.       END IF
  131. *
  132. *     Quick return if possible
  133. *
  134.       IF( N.EQ.0 )
  135.      $   RETURN
  136. *
  137.       IF( N.EQ.1 ) THEN
  138.          IF( ICOMPZ.EQ.2 )
  139.      $      Z( 1, 1 ) = ONE
  140.          RETURN
  141.       END IF
  142. *
  143. *     Determine the unit roundoff and over/underflow thresholds.
  144. *
  145.       EPS = DLAMCH( 'E' )
  146.       EPS2 = EPS**2
  147.       SAFMIN = DLAMCH( 'S' )
  148.       SAFMAX = ONE / SAFMIN
  149.       SSFMAX = SQRT( SAFMAX ) / THREE
  150.       SSFMIN = SQRT( SAFMIN ) / EPS2
  151. *
  152. *     Compute the eigenvalues and eigenvectors of the tridiagonal
  153. *     matrix.
  154. *
  155.       IF( ICOMPZ.EQ.2 )
  156.      $   CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDZ )
  157. *
  158.       NMAXIT = N*MAXIT
  159.       JTOT = 0
  160. *
  161. *     Determine where the matrix splits and choose QL or QR iteration
  162. *     for each block, according to whether top or bottom diagonal
  163. *     element is smaller.
  164. *
  165.       L1 = 1
  166.       NM1 = N - 1
  167. *
  168.    10 CONTINUE
  169.       IF( L1.GT.N )
  170.      $   GO TO 160
  171.       IF( L1.GT.1 )
  172.      $   E( L1-1 ) = ZERO
  173.       IF( L1.LE.NM1 ) THEN
  174.          DO 20 M = L1, NM1
  175.             TST = ABS( E( M ) )
  176.             IF( TST.EQ.ZERO )
  177.      $         GO TO 30
  178.             IF( TST.LE.( SQRT( ABS( D( M ) ) )*SQRT( ABS( D( M+
  179.      $          1 ) ) ) )*EPS ) THEN
  180.                E( M ) = ZERO
  181.                GO TO 30
  182.             END IF
  183.    20    CONTINUE
  184.       END IF
  185.       M = N
  186. *
  187.    30 CONTINUE
  188.       L = L1
  189.       LSV = L
  190.       LEND = M
  191.       LENDSV = LEND
  192.       L1 = M + 1
  193.       IF( LEND.EQ.L )
  194.      $   GO TO 10
  195. *
  196. *     Scale submatrix in rows and columns L to LEND
  197. *
  198.       ANORM = DLANST( 'I', LEND-L+1, D( L ), E( L ) )
  199.       ISCALE = 0
  200.       IF( ANORM.EQ.ZERO )
  201.      $   GO TO 10
  202.       IF( ANORM.GT.SSFMAX ) THEN
  203.          ISCALE = 1
  204.          CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L+1, 1, D( L ), N,
  205.      $                INFO )
  206.          CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L, 1, E( L ), N,
  207.      $                INFO )
  208.       ELSE IF( ANORM.LT.SSFMIN ) THEN
  209.          ISCALE = 2
  210.          CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L+1, 1, D( L ), N,
  211.      $                INFO )
  212.          CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L, 1, E( L ), N,
  213.      $                INFO )
  214.       END IF
  215. *
  216. *     Choose between QL and QR iteration
  217. *
  218.       IF( ABS( D( LEND ) ).LT.ABS( D( L ) ) ) THEN
  219.          LEND = LSV
  220.          L = LENDSV
  221.       END IF
  222. *
  223.       IF( LEND.GT.L ) THEN
  224. *
  225. *        QL Iteration
  226. *
  227. *        Look for small subdiagonal element.
  228. *
  229.    40    CONTINUE
  230.          IF( L.NE.LEND ) THEN
  231.             LENDM1 = LEND - 1
  232.             DO 50 M = L, LENDM1
  233.                TST = ABS( E( M ) )**2
  234.                IF( TST.LE.( EPS2*ABS( D( M ) ) )*ABS( D( M+1 ) )+
  235.      $             SAFMIN )GO TO 60
  236.    50       CONTINUE
  237.          END IF
  238. *
  239.          M = LEND
  240. *
  241.    60    CONTINUE
  242.          IF( M.LT.LEND )
  243.      $      E( M ) = ZERO
  244.          P = D( L )
  245.          IF( M.EQ.L )
  246.      $      GO TO 80
  247. *
  248. *        If remaining matrix is 2-by-2, use DLAE2 or SLAEV2
  249. *        to compute its eigensystem.
  250. *
  251.          IF( M.EQ.L+1 ) THEN
  252.             IF( ICOMPZ.GT.0 ) THEN
  253.                CALL DLAEV2( D( L ), E( L ), D( L+1 ), RT1, RT2, C, S )
  254.                WORK( L ) = C
  255.                WORK( N-1+L ) = S
  256.                CALL DLASR( 'R', 'V', 'B', N, 2, WORK( L ),
  257.      $                     WORK( N-1+L ), Z( 1, L ), LDZ )
  258.             ELSE
  259.                CALL DLAE2( D( L ), E( L ), D( L+1 ), RT1, RT2 )
  260.             END IF
  261.             D( L ) = RT1
  262.             D( L+1 ) = RT2
  263.             E( L ) = ZERO
  264.             L = L + 2
  265.             IF( L.LE.LEND )
  266.      $         GO TO 40
  267.             GO TO 140
  268.          END IF
  269. *
  270.          IF( JTOT.EQ.NMAXIT )
  271.      $      GO TO 140
  272.          JTOT = JTOT + 1
  273. *
  274. *        Form shift.
  275. *
  276.          G = ( D( L+1 )-P ) / ( TWO*E( L ) )
  277.          R = DLAPY2( G, ONE )
  278.          G = D( M ) - P + ( E( L ) / ( G+SIGN( R, G ) ) )
  279. *
  280.          S = ONE
  281.          C = ONE
  282.          P = ZERO
  283. *
  284. *        Inner loop
  285. *
  286.          MM1 = M - 1
  287.          DO 70 I = MM1, L, -1
  288.             F = S*E( I )
  289.             B = C*E( I )
  290.             CALL DLARTG( G, F, C, S, R )
  291.             IF( I.NE.M-1 )
  292.      $         E( I+1 ) = R
  293.             G = D( I+1 ) - P
  294.             R = ( D( I )-G )*S + TWO*C*B
  295.             P = S*R
  296.             D( I+1 ) = G + P
  297.             G = C*R - B
  298. *
  299. *           If eigenvectors are desired, then save rotations.
  300. *
  301.             IF( ICOMPZ.GT.0 ) THEN
  302.                WORK( I ) = C
  303.                WORK( N-1+I ) = -S
  304.             END IF
  305. *
  306.    70    CONTINUE
  307. *
  308. *        If eigenvectors are desired, then apply saved rotations.
  309. *
  310.          IF( ICOMPZ.GT.0 ) THEN
  311.             MM = M - L + 1
  312.             CALL DLASR( 'R', 'V', 'B', N, MM, WORK( L ), WORK( N-1+L ),
  313.      $                  Z( 1, L ), LDZ )
  314.          END IF
  315. *
  316.          D( L ) = D( L ) - P
  317.          E( L ) = G
  318.          GO TO 40
  319. *
  320. *        Eigenvalue found.
  321. *
  322.    80    CONTINUE
  323.          D( L ) = P
  324. *
  325.          L = L + 1
  326.          IF( L.LE.LEND )
  327.      $      GO TO 40
  328.          GO TO 140
  329. *
  330.       ELSE
  331. *
  332. *        QR Iteration
  333. *
  334. *        Look for small superdiagonal element.
  335. *
  336.    90    CONTINUE
  337.          IF( L.NE.LEND ) THEN
  338.             LENDP1 = LEND + 1
  339.             DO 100 M = L, LENDP1, -1
  340.                TST = ABS( E( M-1 ) )**2
  341.                IF( TST.LE.( EPS2*ABS( D( M ) ) )*ABS( D( M-1 ) )+
  342.      $             SAFMIN )GO TO 110
  343.   100       CONTINUE
  344.          END IF
  345. *
  346.          M = LEND
  347. *
  348.   110    CONTINUE
  349.          IF( M.GT.LEND )
  350.      $      E( M-1 ) = ZERO
  351.          P = D( L )
  352.          IF( M.EQ.L )
  353.      $      GO TO 130
  354. *
  355. *        If remaining matrix is 2-by-2, use DLAE2 or SLAEV2
  356. *        to compute its eigensystem.
  357. *
  358.          IF( M.EQ.L-1 ) THEN
  359.             IF( ICOMPZ.GT.0 ) THEN
  360.                CALL DLAEV2( D( L-1 ), E( L-1 ), D( L ), RT1, RT2, C, S )
  361.                WORK( M ) = C
  362.                WORK( N-1+M ) = S
  363.                CALL DLASR( 'R', 'V', 'F', N, 2, WORK( M ),
  364.      $                     WORK( N-1+M ), Z( 1, L-1 ), LDZ )
  365.             ELSE
  366.                CALL DLAE2( D( L-1 ), E( L-1 ), D( L ), RT1, RT2 )
  367.             END IF
  368.             D( L-1 ) = RT1
  369.             D( L ) = RT2
  370.             E( L-1 ) = ZERO
  371.             L = L - 2
  372.             IF( L.GE.LEND )
  373.      $         GO TO 90
  374.             GO TO 140
  375.          END IF
  376. *
  377.          IF( JTOT.EQ.NMAXIT )
  378.      $      GO TO 140
  379.          JTOT = JTOT + 1
  380. *
  381. *        Form shift.
  382. *
  383.          G = ( D( L-1 )-P ) / ( TWO*E( L-1 ) )
  384.          R = DLAPY2( G, ONE )
  385.          G = D( M ) - P + ( E( L-1 ) / ( G+SIGN( R, G ) ) )
  386. *
  387.          S = ONE
  388.          C = ONE
  389.          P = ZERO
  390. *
  391. *        Inner loop
  392. *
  393.          LM1 = L - 1
  394.          DO 120 I = M, LM1
  395.             F = S*E( I )
  396.             B = C*E( I )
  397.             CALL DLARTG( G, F, C, S, R )
  398.             IF( I.NE.M )
  399.      $         E( I-1 ) = R
  400.             G = D( I ) - P
  401.             R = ( D( I+1 )-G )*S + TWO*C*B
  402.             P = S*R
  403.             D( I ) = G + P
  404.             G = C*R - B
  405. *
  406. *           If eigenvectors are desired, then save rotations.
  407. *
  408.             IF( ICOMPZ.GT.0 ) THEN
  409.                WORK( I ) = C
  410.                WORK( N-1+I ) = S
  411.             END IF
  412. *
  413.   120    CONTINUE
  414. *
  415. *        If eigenvectors are desired, then apply saved rotations.
  416. *
  417.          IF( ICOMPZ.GT.0 ) THEN
  418.             MM = L - M + 1
  419.             CALL DLASR( 'R', 'V', 'F', N, MM, WORK( M ), WORK( N-1+M ),
  420.      $                  Z( 1, M ), LDZ )
  421.          END IF
  422. *
  423.          D( L ) = D( L ) - P
  424.          E( LM1 ) = G
  425.          GO TO 90
  426. *
  427. *        Eigenvalue found.
  428. *
  429.   130    CONTINUE
  430.          D( L ) = P
  431. *
  432.          L = L - 1
  433.          IF( L.GE.LEND )
  434.      $      GO TO 90
  435.          GO TO 140
  436. *
  437.       END IF
  438. *
  439. *     Undo scaling if necessary
  440. *
  441.   140 CONTINUE
  442.       IF( ISCALE.EQ.1 ) THEN
  443.          CALL DLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV+1, 1,
  444.      $                D( LSV ), N, INFO )
  445.          CALL DLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV, 1, E( LSV ),
  446.      $                N, INFO )
  447.       ELSE IF( ISCALE.EQ.2 ) THEN
  448.          CALL DLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV+1, 1,
  449.      $                D( LSV ), N, INFO )
  450.          CALL DLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV, 1, E( LSV ),
  451.      $                N, INFO )
  452.       END IF
  453. *
  454. *     Check for no convergence to an eigenvalue after a total
  455. *     of N*MAXIT iterations.
  456. *
  457.       IF( JTOT.LT.NMAXIT )
  458.      $   GO TO 10
  459.       DO 150 I = 1, N - 1
  460.          IF( E( I ).NE.ZERO )
  461.      $      INFO = INFO + 1
  462.   150 CONTINUE
  463.       GO TO 190
  464. *
  465. *     Order eigenvalues and eigenvectors.
  466. *
  467.   160 CONTINUE
  468.       IF( ICOMPZ.EQ.0 ) THEN
  469. *
  470. *        Use Quick Sort
  471. *
  472.          CALL DLASRT( 'I', N, D, INFO )
  473. *
  474.       ELSE
  475. *
  476. *        Use Selection Sort to minimize swaps of eigenvectors
  477. *
  478.          DO 180 II = 2, N
  479.             I = II - 1
  480.             K = I
  481.             P = D( I )
  482.             DO 170 J = II, N
  483.                IF( D( J ).LT.P ) THEN
  484.                   K = J
  485.                   P = D( J )
  486.                END IF
  487.   170       CONTINUE
  488.             IF( K.NE.I ) THEN
  489.                D( K ) = D( I )
  490.                D( I ) = P
  491.                CALL DSWAP( N, Z( 1, I ), 1, Z( 1, K ), 1 )
  492.             END IF
  493.   180    CONTINUE
  494.       END IF
  495. *
  496.   190 CONTINUE
  497.       RETURN
  498. *
  499. *     End of DSTEQR
  500. *
  501.       END
  502.